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I. INTRODUCTION 

Properties of geometrically frustrated spin systems in 
various dimensions, geometries, and temperature regimes 
are at the heart of modern condensed matter physics. 
Here, frustration is a technical term which refers to the 
presence of competing forces that cannot be simultane- 
ously satisfied. In numerous quantum antiferromagnets 
frustration often has a simple geometric origin. Localized 
spins on two- and three-dimensional lattices with triangu- 
lar motifs, such as planar triangular antiferromagnet and 
three-dimensional pyrochlorc antiferromagnets, cannot 
assume energetically favorable antiparallel orientation. 
Of three spins forming a minimal triangle, and interact- 
ing via simple antiferromagnetic pair-wise exchange in- 
teraction, only two can be made antiparallel, leaving the 
third one frustrated. In case of classical Ising spins, which 
can point up or down with respect to some axis, this leads 
to an extensive ground-state degeneracy: for example, in 
a system of N Ising spins on a triangular lattice there 
are Ojv = e°' 323Ar configurations having the same (min- 
imal) energy. On a three-dimensional pyrochlore lattice 
of site-sharing tetrahedra, the (effectively) Ising spins of 
Dy 2 Ti 2 07, H02TI2O7 and Ho 2 Sn 2 07 realize 1 fascinat- 
ing spin ice physics 2 where strong local ice rules (for any 
given tetrahedron, two of its spins must point in, and the 
other two - out) enforce long-ranged power-law correla- 
tions between spins 3 , in effect realizing artificial magnetic 
field and fractionally charged magnetic monopolcs 4 ! 

Quantum spins can exploit this extensive degeneracy 
via quantum-mechanical coupling between different con- 
figurations - their wave function can be thought of as a 
linear superposition of all degenerate microstates repre- 
sented by classical patterns of up- and down-spins. In the 
case of strong coupling we may arrive at a quantum spin 
liquid (QSL) 5 ' 6 state in which spins never settle in one 
particular configuration and continue their exploration 
forever. It is clear that such a state encodes highly non- 
trivial correlations between different spins when flipping 



of one spin induces that of its neighbors so that as a 
whole the spin system remains in the lowest-energy man- 
ifold. Extensive experimental 7-20 and theoretical 21-30 
search for materials and models which may realize this 
intriguing QSL state constitutes one of the main topics 
of the quantum spin physics. Currently there are several 
intensely researched materials that hold promise of real- 
izing the elusive spin liquid state. Among them, we men- 
tion two-dimensional spin- 1/2 organic triangular antifer- 
romagnets EtMe 3 Sb[Pd(dmit)2]2 10,11 ' 19 and k-(BEDT- 
TTF) 2 Cu 2 (CN) 3 8 ' 19 , spin-1 material NiGa 2 S 4 16 , a se- 
ries of inorganic quasi-two-dimensional kagome lattice 
antiferromagnets: herbersmithite ZnCu3(0H)6Cl2 12 ' 15 , 
volborthite Cu 3 V 2 7 (OH) 2 -2H 2 7 ' 13 , and vesignieite 
BaCu3V 2 0g(0H) 2 14 ' 18 , and a three-dimensional hyperk- 
agome antiferromagnet Na4ir30s 9 . 

Taking the system of frustrated spins to a finite tem- 
perature, where all experiments are done, adds thermal 
randomness to the picture. At a finite temperature T 
even classical spins can explore different microstates from 
the lowest-energy manifold. This too leads to a strongly 
correlated (although not necessarily phase-coherent, as 
in the case of quantum spins at T = 0) motion of spins 
which is often described by a term "cooperative param- 
agnet" . Even if the ground state of spin system is not a 
true spin liquid, but instead is one of the many possible 
ordered states, the spins will (thermally) disorder at suf- 
ficiently high temperature, T > T , where To stands for 
the ordering temperature. In the usual, non-frustrated 
magnets the ordering temperature is determined by the 
exchange interaction energy J and coordination number 
of the lattice z, T ~ zS(S + 1) J. Quite generally, it 
is of the order of the Curie- Weiss temperature CW which 
is easily determined experimentally via high-temperature 
behavior of the spin susceptibility \ cx (T — ^ cw ) -1 . In 
antiferromagnets 9 CW is negative, and, in the absence of 
frustration, its absolute value sets the scale at which 
correlations between spins become pronounced. Thus, 
Tq ~ |0 C w|- Frustrated magnets are very different as 
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there Tq -C |0 CW | : despite experiencing strong interac- 
tions with each other, the spins can not "agree" on one 
particular pattern which would satisfy them all. In the 
case of a true spin liquid Tq = and the order never 
arrives. In the majority of studied situations the order 
does take place, To > 0, but only at a temperature much 
lower than the naive estimate provided by \0 cvr \. This 
results in a wide temperature interval Tq < T < \9 CW \ 
where the spins are strongly interacting but remain in 
a disordered cooperative paramagnet state. In fact, the 
existence of such temperature (and energy) window rep- 
resents a defining feature of the frustrated magnet, as 
argued by Ramirez 31 who introduced the frustration pa- 
rameter / = |# cw |/To (it is not uncommon to find a sit- 
uation with / ~ 100 or greater). 

Strong suppression of the ordering temperature can be 
also the consequence of close proximity to the quantum 
critical point, separating the spin-disordered state (such 
as a putative QSL) from the more usual ordered one. 
It is well established now 32 that the finite-temperature 
region above the quantum critical point, known as the 
quantum- critical region, is as informative of the quantum 
state of the many-body system as the unreachable T = 
ground state. The quantum- critical scaling of, say, the 
dynamic spin susceptibility contains information on the 
spin correlation length, dynamic exponent z and other 
critical exponents characterizing spin systems of different 
symmetries and dimensions. 

Importantly, there is a large number of high-precision 
probes — neutron and X-ray scattering, nuclear magnetic 
resonance, muon spin rotation, susceptibility, magneti- 
zation and specific heat measurements — which allow us 
to address various aspects of strange and conflicting be- 
havior of frustrated quantum magnets experimentally in 
a wide range of energies and temperatures. Given that 
in many important cases the correlated spin-liquid re- 
gion Tq < T < \9 CW \ occupies most of the experimentally 
accessible temperature interval, unbiased understanding 
of correlations and dynamics in this regime becomes the 
major theoretical task. 

Theoretical understanding of frustrated magnetism at 
finite temperature is severely limited by the lack of nat- 
ural small parameter (s). As a result, possible analyti- 
cal approaches require one to study suitably 'deformed' 
models, such as, for example, very popular and well de- 
veloped large-component (large- N) version of the Heisen- 
berg spin model on frustrated lattices. The small pa- 
rameter is then provided by 1/N, expansion in powers 
of which (about the N = oo limit) controls the calcu- 
lation. However the physical limit of S — 1/2 SU(2) 
lattice spins corresponds to rather small N: 3 in the 
case of O(N) generalization and 1 in the case of Sp(N). 
Whether or not continuation of the results from N = oo 
to the physical value is reliable remains an open (and 
case sensitive) question. Other popular 'deformations' 
include SU(2)-to-U(l) symmetry reduction(s) and quan- 
tum dimer model approaches. While very insightful and 
interesting in their own, applicability of these 'modifica- 



tions' to the original problem is always an issue. While 
analytical approaches are extremely useful in providing 
us with qualitative physical insights and understanding, 
in frustrated magnets they often fall short of quantitative 
description desired by experimentalists. Numerical ap- 
proaches are limited as well. Exact diagonalizations are 
restricted to small systems (about 40 sites at most) due 
to the exponentially large Hilbert space while standard 
quantum Monte Carlo techniques usually suffer from the 
infamous 'sign problem'. Powerful series expansion meth- 
ods often start to diverge in the most interesting regime 
Tq < T < \0 CW \. Variational tensor-network type meth- 
ods are also suffering from finite-size limitations and are 
mostly limited to ground state properties. 

In this article we combine the most versatile theoretical 
tool, Feynman diagrammatics, with the power of Monte 
Carlo sampling of complex configuration spaces. The 
simplest way of arriving at the diagrammatic technique 
for spins is to represent them by auxiliary fcrmions with 
imaginary chemical potential. This trick was introduced 
by Popov and Fedotov for spin- 1/2 systems in Ref. 33. 
The ultimate strength of the diagrammatic approach as 
compared, e.g. with high-temperature or strong coupling 
expansions, comes from its self-consistent (skeleton) for- 
mulation with automatic summation of certain classes 
of graphs up to infinite order. This lead to better con- 
vergence properties and the possibility of obtaining re- 
liable results in the strong coupling regime. It turns 
out that skeleton formulations can be easily implemented 
within the sampling protocols leading to the so-called 
Bold Diagrammatic Monte Carlo (BDMC) 34 which ob- 
tains physical answers by computing contributions from 
millions of graphs and extrapolates them to the infinite 
diagram order. Recently, this method was successfully 
applied to the normal state of the Fermi-Hubbard model 
at moderate interaction strength 35 and the strongly cor- 
related system of unitary fermions (the so-called BCS- 
BEC crossover problem) 36 . It was, however, never im- 
plemented for models of quantum magnetism, which, at 
least at the formal level, can be also viewed as a fermionic 
system with strong correlations. Here we present the first 
attempt to achieve an accurate theoretical description of 
the correlated paramagnetic regime within the BDMC 
framework. 

In what follows in Sec. II we consider a spin- 1/2 Heisen- 
berg type model at finite temperature and its auxiliary 
fermion version which admits the standard diagrammatic 
expansion. In Sec. Ill we formulate principles of the 
BDMC technique and the specific self-consistent scheme 
for dealing with skeleton diagrams based on fully dressed 
lines. We proceed with detailed description of the worm 
algorithm for efficient sampling of the resulting config- 
uration space (updates, counters, and data processing) 
in Sec. IV. Our results for the triangular lattice anti- 
ferromagnet are presented, discussed, and compared to 
the best available finite temperature simulations based 
on numerical linked-cluster (NLC) expansions 37 and ex- 
trapolations of high-temperature series 38 in Sec. V. We 
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conclude with broader implications of this work and fur- 
ther developments in Sec. VI. 



II. MODEL AND ITS 'BOLD-LINE' 
DIAGRAMMATIC EXPANSION 

Consider the standard Heisenberg model 

H = 22 Jij Si ■ Sj , (1) 

i,j 

where S are quantum spin-1/2 operators. The dimension 
of space, lattice geometry, and interaction range are as- 
sumed to be arbitrary. Now, the idea is to replace spin 
degrees of freedom with fermionic ones: 



(2) 



a,/3 



Here fip is the second quantized operator annihilating 
a fermion with spin projection j3 = ±1 on site i, and 
<7 are the Pauli matrices. As a result, we convert cou- 
pling between spins into the standard two-body inter- 
action term J ijta p lS fLfipfjyfjs witn matrix element 
Jij,af3y8 = (1/4) Jij&a0 • &-)8- By doing so we also in- 
crease the Hilbert space on every site from 2 to 4 by 
adding non-physical states with zero and two fermions. 

The benefit of the auxiliary fermion representation 
does not require explanation: once the spin model is 
mapped onto a familiar problem of interacting fermions 
one can employ numerous diagrammatic tricks to solve it. 
However, this raises an issue of eliminating contributions 
from unphysical states to the answer in a manner con- 
sistent with the diagrammatic technique. Remarkably, 
there is a very simple way to achieve the goal by adding 
the chemical potential term to the fermionic Hamilto- 
nian 33 (see also Ref. 39) 



H F = 



Jij,a/3yS fiafipfjryfjg 



1), (3) 



with complex /j, = -inT/2 and n, ; = J2 a fLfta- Tnc 
added term commutes with the original Hamiltonian 
and has no effect on properties of the physical subspace 
{n, = 1} whatsoever. Moreover, the grand canonical 
partition functions and spin-spin correlation functions of 
the original spin model and its fermionic version are also 
identical because (i) physical and non-physical sites de- 
couple in the trace and (ii) the trace over non-physical 
states yields identical zero on every site. As a result, we 
arrive at a rather standard Hamiltonian for fermions in- 
teracting through two-body terms. Complex value of the 
chemical potential is essentially a zero price to pay for 
the luxury of having the diagrammatic technique. 

Formally, the entire setup is similar to the fractional 
quantum Hall effect system because the non-interacting 



part of Eq. (3) describes particles with zero dispersion re- 
lation, i.e. we start building the solution from the degen- 
erate manifold. There is one important difference though: 
in conventional fermionic systems both the Green's func- 
tion and the spin-spin (or density-density) correlation 
function contain direct physical information about the 
system. In spin systems the Green's function is rather 
an auxiliary object which always remains localized on a 
single site. This does not imply any pathological behav- 
ior yet since the physical degrees of freedom are spins, 
i.e. bilinear combinations of fermionic operators. Corre- 
spondingly, the main object of interest is not the Green's 
function of the system but the spin-spin correlation func- 
tion, or magnetic susceptibility: 

X (i,j,r) = (T T S?(0)SJ(r)) = ±{T T §i(0) ■ S^r)) , (4) 

where T T stands for the imaginary time ordering opera- 
tor. 

To simplify the presentation below we consider lattices 
with one atom per unit cell and make use of the lattice 
translation invariance; otherwise one would need to keep 
an index enumerating different sites in the unit cell. For 
the same reason we do not place the system into the ex- 
ternal magnetic field to preserve the symmetry between 
up- and down-spins; in the presence of the magnetic field 
Hq one has to add spin-dependent real part to the chem- 
ical potential /i — ¥ fi a = —iirT/2 — aHo and take proper 
care of the spin index. 

The diagrammatic technique itself for Eq. (3) is ab- 
solutely standard 40 . The perturbative diagrams are 
expressed in terms of the non-interacting, or 'bare', 
Green's functions (particle propagators), G^l(i, r) = 

SiflS a pG^(r), and two-body interaction lines Ji-j^p-yS- 
The non-interacting Green's function on a lattice in the 
imaginary time representation reads (below the Boltz- 
mann constant kg = 1) 



G (0) (t >0) = - 



l + cxp^i/T) i- 1 ' 



(5) 



with conventional anti-periodic boundary conditions 
G(°)(t < 0) = -G(°)(l/T + r). In what follows we com- 
pletely suppress the site index for purely local quantities. 
For the Heisenberg model the dependence of the interac- 
tion line on spin indexes is rather simple 



J, 



Ji 



1-3 



M, 



M a pjs = aj 5 a ,p8~f t s + 2 8s t0l 8^- a 5p- a , (6) 

The first term describes diagonal coupling between the 
spin densities on sites i and j, while the second spin- flip 
term exchanges spin values, see Fig. 1. The magnitude 
of the spin-flip process is fixed by the SU(2) symmetry 
of the problem. Even when interactions arc screened by 
many-body effects, see below, the retarded interaction 
has exactly the same dependence on spin indexes. An 
arbitrary perturbative diagram of order n contributing 
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to, say, the free-energy is obtained by (i) placing graph- 
ical elements depicted in Fig. 1 with some space/time 
variables i, j, and r € (0, and connecting incoming 
and outgoing propagators with the same spin and site 
index to each other in such a way that all points in the 
resulting graph are connected by some path. 




Diagonal 



-J tJ S(T 2 -Tj 



Spin-flip 



FIG. 1. (Color online) Graphic representation of allowed 
interaction processes for the Heisenberg model. 



III. BOLD DIAGRAMMATIC MONTE CARLO 
SCHEME 

The unique feature of diagrammatic expansions for 
propagators is that there are no numerical coefficients 
in the diagram weight depending on the diagram order 
or structure (this is not true for other well-known series 
such as virial, high-temperature, linked-cluster, strong- 
coupling, etc. expansions). This leads to the diagram- 
matic technique when certain infinite sets of diagrams, 
e.g. in the form of geometric series, are easily dealt with 
by algebraic means or reduced to self-consistently defined 
integral equations. 

In the skeleton technique, the diagrams are classified 
according to some rule which eliminates the need for 
computing repeated blocks of diagrams. In the simplest 
scheme which is used in this paper, one identifies the 
proper self-energy blocks which consist of diagrams where 
all vertexes (points where two particle propagators and 
the interaction line meet) remain connected by some path 
when one removes any two lines of the same kind: two 
propagator lines with the same spin index or two inter- 
action lines. We will refer to this set of diagrams as 
irreducible. The omitted diagrams are fully accounted 
for by replacing bare propagators and interaction lines 
in the proper self-energy blocks with exact propagators, 
G(t), and screened interactions, W(r, t)M . The result- 
ing formulation is self-consistent and highly non-linear 
since G and W depend on proper self-energies through 
the Dyson type equations. If £ is the proper self-energy 
for the particle propagator, and IT is an analogous quan- 
tity for the interaction line (better known as polarization 
operator) then (in Fourier representation (r, r) — > (q,uj m ) 



for space-time variables) 
G(°)(m) 



G(m) 



W(q, m) 



G(°)(m)£(m) 



W(q,m). (7) 



4 - J(q)ll(q, m) 4 

where J(q) = J2 r e zqr J(r). Due to fermionic/bosonic na- 
ture of propagators G/W we have different definitions 
of the Matsubara frequency here, w m = 2-KT(m + 1/2) 
for (G,£) and w m = 2nTm for (W,II). Note also that 
we split the W function into two parts by separating 
out the original coupling. This is done for technical rea- 
sons explained below; here we simply point out that in 
the imaginary time domain this is equivalent to pay- 
ing special attention to the ^-functional contribution, 
W(q,r) — 4 J(q) 6(t) + W(q,r). Correspondingly, for 
graphical representation of the diagrams we use wavy 
lines for W and dashed vertical lines for bare coupling. 
The corresponding (rather standard in many-body the- 
ory) G 2 VF-skeleton scheme is illustrated in Fig. 2. The 
magnetic susceptibility Eq. (4) is directly related to the 
polarization operator appearing in Eq. (7) and Fig. 2 

Il(q,n) 



4- J(q)n(q,i 



(8) 



w 





a/3 



FIG. 2. (Color online) Typical low-order diagrams contribut- 
ing to the particle self-energy and polarization operator within 
the G 2 W / -skeleton scheme (red color denotes spin-exchange 
coupling). Diagrams for E and IT, in their turn, are used to 
calculate fully dressed G and W functions, see Eq. (7). 

To make connection with general rules of diagrammatic 
MC 41 ' 42 we formalize the problem at hand as computing 
quantity Q(y, s) (where y stands collectively for space, 
imaginary time, and spin variables, while s = 1,2 labels 
the proper self-energy and the polarization operator sec- 
tors) from the series of multidimensional sums/integrals 

Q(y,s) = J2 f d Xl ...dx n dY V(n^,{x i };Y,s)S(y -Y) . 

(9) 
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Here n = 0, 1, ... cm is the diagram order, £ labels dif- 
ferent terms of the same order, internal integra- 
tion/summation variables, and T> is the diagram contri- 
bution to the answer. Formally, one can think of the 
set of skeleton diagrams for the free-energy of the system 
with one line marked as 'dummy'. When the dummy line 
is removed from the graph the rest is interpreted as a di- 
agram for Ti(y) (in this case s = 1 and the dummy line 
is the particle propagator) or n(y) (in this case s = 2 
and the dummy line is the W one). The last rule fol- 
lows from the fact that II(y) is a continuous function 
of time. In the BDMC approach one interprets Eq. (9) 
as averaging of e iaTS ^ D ' S(y — Y) over the configuration 
space v = (n, £, Xi, . . . , x n ; Y, s) with probability density 
proportional to \T>\. 

The skeleton formulation does not cause any funda- 
mental problem for Monte Carlo methods and is easy to 
implement for diagrams of arbitrary order. Essentially, 
at any stage in the calculation both G and W are con- 
sidered to be known functions (the calculation may start 
with G = G(°) and W = J /A) while Eq. (7) is used from 
time to time to improve one's knowledge about G and W 
using accumulated statistics for E and II and fast Fourier 
transform algorithms. Moreover, we have shown 34 that 
BDMC methods are more stable and have better con- 
vergence properties than conventional iterations. One 
immediately recognizes that in the skeleton-type formu- 
lation (i) the number of diagrams to be sampled in a given 
order is dramatically reduced, (ii) the convergence of the 
skeleton series is likely to be different than that of the 
bare series, (iii) non-analytic and non-perturbative be- 
havior might emerge even from a finite number of terms 
due to highly non-linear self-consistent formulation (we 
refer here to the famous mean-field BCS solution). 

The other crucial advantage of BDMC over more con- 
ventional MC methods simulating finite clusters of spins 
is that it deals directly with the thermodynamic limit 
of the system. In practice, in dimension d > 1 the 
error bars often become too large before a reliable ex- 
trapolation to the thermodynamic limit can be done. In 
this sense, BDMC is not subject to the infamous sign- 
problem which is understood as exponential scaling of 
computational complexity with the space-time volume of 
the physical system. Feynman diagrams do alternate in 
sign and contributions from high-order diagrams cancel 
each other to near zero. However, this behavior is better 
characterized as a 'sign-blessing', not a sign-problem, be- 
cause it is crucial for convergence properties. With the 
number of graphs growing factorially with their order the 
only possibility for obtaining series with finite conver- 
gence radius is to have sign-alternating terms such that 
high-order diagrams cancel each other (the sign-blessing 
phenomenon). For series with finite convergence radius 
there are numerous unbiased re-summation techniques 
which allow one to determine the answer well outside of 
the convergence radius provided enough terms in the se- 
ries are known. Relatively small configuration space for 
skeleton diagrams allows one to establish if sign-blessing 



takes place in a given model and to obtain accurate re- 
sults for diagram orders as large as 7-10 (depending on 
the model). 

It should be noted that in recent years the Popov- 
Fedotov trick has become popular within the functional 
renormalization group (PFFRG) framework. It has 
been applied to frustrated J\ — J2 model on square 
lattice 43 , planar J\ — J 2 — J3 antiferromagnet 44 , spa- 
tially anisotropic triangular antiferromagnet 45 , and hon- 
eycomb lattice antiferromagnet with competing interac- 
tions 46 . These studies also attempt at attacking the 
problem using Feynman diagrammatic series but are rad- 
ically different in the technical implementation. While 
PFFRG is also based on sums of subsets of selected di- 
agrams to infinite order, it does not offer a convenient 
way to check for convergence of final results as more 
and more diagrams are retained. This is the most severe 
drawback of PFFRG; after all the major problem with 
existing theories is reliable estimate of the accuracy. The 
PFFRG spectral functions are very broad in energy 43 
and appear to underestimate ordering fluctuations. This 
also leads to significant rounding of susceptibility peaks 
and makes identification of different phases difficult. In 
addition, these studies are typically focused on the zero- 
temperature phase diagram, not the finite-temperature 
cooperative paramagnet state. 

IV. NORMALIZATION AND 
WORM-ALGORITHM UPDATES 

To simplify notations let us write the diagrammatic 
contribution from the configuration space point v as 
T> 1, = e %Vu D v and call the non- negative function D v the 
configuration 'weight'. Within the G 2 VF-skeleton scheme 
D u is given by the modulus of a product which runs over 
all lines 

D v = I J] /line (") I , (10) 
lines 

where f\i ne (v) stands for a collection of functions describ- 
ing various lines in the diagram. At this point we notice 
that equation 

Q(y,s) = ^2e^D v 6(y-Y) . (11) 

can be always interpreted as averaging over the proba- 
bility density distribution P vs = D U /C S1 where C s is the 
normalization factor, and thus sampled by MC methods: 

Q(y, s)=C s Y / e^6(y-Y)P vs — ► C s £ e*"5(y-Y) . 

(12) 

In the last transformation we replace the full sum over 
the configuration space with the stochastic sum over con- 
figurations which are generated from the probability den- 
sity P vs . This is, of course, nothing but the standard MC 
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approach to deal with complex multi-dimensional spaces. 
We stress here that all configuration parameters are sam- 
pled stochastically, including the diagram order and its 
structure, making Diagrammatic MC radically different 
from calculations which first create a list of all diagrams 
up to some high-order and then evaluate them one-by- 
one (often with the use of MC methods for doing the 
integrals). 



magnetic field the spin-up and spin-down contribution 
exactly cancel each other. Correspondingly, 

MC 

■Zl = <5^,(Hartrcc,s=l) i 

v 
MC 

^2 = ^2$u,(n=l,s=2) ■ (16) 



Normalization 



B. Worm-algorithm trick 



The normalization constant C can be determined in a 
number of ways: 

(i) using known behavior of Q(y, s) in some limiting case, 
for example Q(y y , s) -> Qo(s), 

(ii) through the exact sum rule, J dyQ(y,s) — R s , if 
available, or, more generically, 

(iii) by measuring the ratio between the contributions of 
all diagrams in Eq. (9) and diagrams which are known 
either analytically or numerically with high accuracy. In- 
deed, imagine that one or several diagrams, say the lowest 
order ones, are known and their integrated contribution 
to the answer is Qn(s) = J dyQw(y, s). Let J\f s be their 
configuration space. Then the ratio Q{y,s)/Q^{s) can 
be measured in the MC simulation as 



dyQ(y,s) 
This leads to 



MC 



MC 



/Q N (s) = £>^)/(£<W s e^) . 

(13) 



: N 



(s) 



MC 



(14) 



The diagrams used for normalization are not necessar- 
ily the physical ones, i.e. they can be artificially "de- 
signed" to have simple analytic structure and added as a 
special sector to the configuration space {i/}, see Ref. 34. 
In the latter case, the diagrams contributing to Q(y, s) 
and Z s are mutually exclusive and physical contributions 
have to be filtered by (1 — S^^J. 

In the present study we use the modulus of the Hartree 
diagram to normalize statistics for S, and the modulus of 
the lowest order GG-diagram, see the first term in Fig. 2, 
to normalize statistics for II: 

g N (l) = E Ar = ^^|^||G(r = -0)|, 

r a 

Q N (2) = n N = V [ 1/T dr\G(T)G(-r)\ . (15) 

Even though in the self-consistent scheme one does not 
know the G(r)-function analytically (it is tabulated nu- 
merically) it takes no time to compute the normalization 
factor Iljv with high accuracy. We take the modulus of 
the Hartree diagrams because in the absence of external 



We now proceed with the description of updates which 
ensure that all points in the configuration space are sam- 
pled from the probability density P vs . Among many pos- 
sibilities we seek a scheme which involves the smallest 
number of lines in a single update and does not require 
global analysis of the diagram structure. Such updates 
are called "local". Typically, they are much more flexi- 
ble in design, are easier to implement, and lead to more 
efficient codes by having large acceptance ratios. The 
reader not interested in algorithmic details may proceed 
directly to the next Section. 

An easy way to ensure that the diagram is irreducible 
is to check that no two lines in the graph have the same 
momentum — by momentum conservation laws the iso- 
lated self-energy blocks cannot change the line momen- 
tum. By creating a hash table where momenta of all lines 
are registered according to their values one can readily 
verify that momenta of updated lines are not repeated in 
the graph without looking at the graph topology or ad- 
dressing all other lines. This simple tool solves the prob- 
lem of performing local updates within the irreducible 
set. It can be applied even if diagrams are sampled in 
the real-space representation, as is done in this article, by 
attaching an auxiliary momentum variable to each line 
and satisfying the momentum conservation law at each 
vertex. In this case, the sole purpose of introducing aux- 
iliary momenta is an efficient monitoring of the diagram 
topology and the T> v value is independent of them. 

Since the trick with auxiliary momenta is based on mo- 
mentum conservation laws one is necessarily limited to ei- 
ther (i) performing updates on closed loops, i.e. partially 
abandoning the idea of local updates, or (ii) extending 
the configuration space to include diagrams which vio- 
late these conservation laws. The second strategy is the 
essence of the worm algorithm. One more reason for us- 
ing the worm algorithm approach is the spin projection 
conservation law in the interaction process, see Eq. (6) 
and Fig. 1. It also requires that updates are performed 
only on closed loops of interaction lines and propagators, 
unless one admits diagrams which violate the correspond- 
ing conservation law. It turns out that a straightforward 
extension of the worm algorithm introduced in Refs. 47 
and 35 allows one to go around both hurdles. 

The additional (unphysical) diagrams have the follow- 
ing structure. There are two special vertexes, or 'worms', 
S and T, where momentum and spin conservation laws 
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are violated, see Fig. 3. In a given graph, any two ver- 
texes can be special if they are not connected by the in- 
teraction line. Conservation laws become satisfied if one 
imagines a special line connecting S to T. This line is 
not associated with any propagator and its mission is to 
transfer momentum p w and spin projection 1 (doted line 
in Fig. 3); otherwise special vertexes represent a drain 
and a source of momentum p w and spin projection 1, see 
Figs. 3 and 4. It has to be realized that the spin conser- 
vation law is formulated for a pair of vertexes connected 
by the interaction line, i.e. an interaction line contain- 
ing a worm on one of its ends is unphysical and can be 
described by any suitable function F(r,r). 



/ S 



(P„..* : =D 




FIG. 3. (Color online) Non-physical diagram with two 
special vertexes (worms) 5 and T (marked by blue circles). 
Momentum and spin conservation laws would be satisfied at 
all vertexes and interaction lines if one considers S/T as a 
drain/source of momentum p w and spin projection 1; the same 
rule is recovered by imagining a line (green dotted) which car- 
ries (p m ,s z = 1) from 5 to T. If special vertexes were not 
present in the diagram, then this graph would contribute to 
E(ti — T2) after removing the dummy propagator line marked 
by the cross (if cross marks the dummy W^-line then the dia- 
gram is contributing to II). 

The worm algorithm idea is based on the observation 
that updates performed with the use of special vertexes 
can always be made local, including non-trivial changes 
in the diagram topology and order as well as transforma- 
tions replacing diagonal interaction lines with spin-flip 
ones. Even though unphysical diagrams are frequently 
encountered in the simulation process they are excluded 
from the statistics of E and II which is accumulated only 
on the physical set of diagrams. 



C. Updates 

Below we describe the simplest ergodic updating 
scheme. With trivial modifications and additional filters 
for proposals which would be rejected because they are 
incompatible with the allowed configuration space it can 
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FIG. 4. (Color online) Detailed structure of S and T ver- 
texes. Note, that the interaction line containing a worm on 
one of its ends is also unphysical and the spin conservation 
law applies to a pair of vertexes. 



be made more efficient. This, however, would overwhelm 
the presentation with minor programming details and we 
choose not to discuss them here. For example, below we 
will use the following algorithmic rules for dummy lines 
used to identify diagrams as S- or II-type: (i) the dummy 
interaction line cannot be removed or created in any up- 
date, (ii) if the dummy propagator originating from ver- 
tex A is modified by adding/removing an intermediate 
vertex C then A always remains the originating vertex 
of the new dummy line. These rules can be easily mod- 
ified to avoid fast rejections of updates in certain cases 
at the expense of using additional random numbers to 
deal with available choices and taking care of them in ac- 
ceptance ratios. Alternatively, the notion of the dummy 
line can be avoided altogether by designing an improved 
estimator based on free-energy diagrams. 

Create - Delete The pair of complementary updates 
Create - Delete switches between physical and unphysical 
sectors by inserting/removing a pair of special vertexes 
connected by the particle propagator. It is not allowed 
to have S = T or to have S and T being connected 
by the interaction line, see an illustration in Fig. 5. In 
Create, the particle propagator for update and the type 
of special vertex to be placed at its origin [A S or 
A — > T) are selected at random. One has to verify that 
flipping the propagator spin is consistent with the worm 
rules or reject the proposal. The missing momentum p w 
at the worm vertex is selected at random. In Delete, one 
selects a worm at random, checks that the outgoing par- 
ticle propagator arrives at the other worm, and proposes 
to remove worms from the diagram. The corresponding 
acceptance ratios are given by 

-Rcroatc = ~p~ 2n , i?Dclctc = ~ — , (17) 

U u u\ D v 1*2 In 

where n is the diagram order and v and v' are configu- 
ration space points before and after the update and D v 
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is the product in Eq. (10). In the present case only one 
propagator and two interaction lines are affected: 



D„ 



GAB{ot)FscFTD 



GABi-aWAcWi 



BD 



(18) 



in Create and similarly, D v >/D v = \GWW/GFF\, in 
Delete with appropriate arguments for all functions in- 
volved. [If the line is labeled as W it can be either of 
J or W type.] In what follows we will stop mentioning 
which vertexes determine function parameters since these 
can be easily recovered from the figures. Finally, assum- 
ing that the protocol of deciding which update should 
be implemented next is random and based on assigning 
each update some probability, u^, we mention the ratio 
of probabilities in the acceptance ratio (u\ for Create and 
U2 for Delete). 
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FIG. 5. Two cases for Create and Delete updates which in- 
sert/remove a pair of worms at the ends of the particle prop- 
agator. 

Create-H - Delete-H This pair of complementary up- 
dates also switches between physical and unphysical sec- 
tors with an additional ingredient — it increases/decreases 
the diagram order by attaching a Hartree-type bubble to 
the existing graph. According to the rule 'no two lines 
may have the same momentum' the diagram remains irre- 
ducible because one of the worms is placed on the bubble 
vertex. The overall transformation is illustrated in Fig. 6; 
the text below addresses to this figure with regards to 
the procedure of selecting specific graph parameters. In 
Create-H a particle propagator (going from vertex A to 
vertex B) and whether to place S or T on vertex B is de- 
cided at random. If the proposal is inconsistent with the 
worm rules it has to be rejected. Next, a new time vari- 
able for the intermediate vertex C is generated from the 
probability density, t(r), and a random decision is made 
whether the new interaction line is of the J- or IF-type. 
For the J-line (with r' = r) the position of the second 
worm vertex in space is obtained from the normalized 
X(r') oc J(r') distribution. For the IF-line this position 
is obtained from some designed probability distribution 



Y(r') while the time location is drawn from the prob- 
ability density, i(r'). The spin variable in the bubble, 
the bubble momentum variable p and the worm momen- 
tum p w are decided at random. In Delete-H a random 
choice is made what type of special vertex must be on 
the Hartrec bubble provided the overall topology of lines 
is identical to that on the r.h.s of Fig. 6. The proposal 
is to remove worms and the bubble from the diagram. It 
is rejected if either the propagator originating from C or 
the interaction line attached to C is the dummy one. The 
acceptance ratios for these updates are (probabilities of 
calling Create-H are Delete-H are u 3 and W4, respectively) 



Rc 



rcatc— H 



D v u 3 t(r) 



l/X(r'); 
l/Y(r')t(r'); 



(J) 
(W) 



(19) 



clctc— H 



IV u 3 t(r) 
D v u 4 2 3 (n - 1) 



X(r<); (J) 
Y(r')t(r'); (W) 

with the diagram weight ratios D v i / D v given by 
\GGFFG/GW\ and \GW/GGFFG\ in Create-H and 
Delete-H, respectively (note that here n is the initial di- 
agram order). The simplest choices for probability dis- 
tributions in Eqs.(19) and (20) would be uniform distri- 
butions t(r) = T, X(r') = 1/z, and Y(r') = 1/V, where 
V is the total number of lattice sites. One can use other 
functions for better acceptance ratio. The three updates 
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FIG. 6. Two cases for Create-H and Delete-H updates which 
insert / remove a pair of worms and increase the diagram order 
by adding a Hartree-type bubble. 

described next represent a random diffusion of special 
vertexes along the graph lines (this places S and T on 
any allowed pair of vertexes) supplemented by an update 
which changes the graph topology. 

Move-P In this self-complementary update one se- 
lects at random S or T and proposes to shift the se- 
lected worm along the incoming or outgoing propagator 
line, deciding again randomly between the two choices. 
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Since all four case are identical in their implementation 
we describe below an update shifting S along the outgo- 
ing propagator line to vertex B, see upper panel in Fig. 7, 
if B 7^ S, of course. According to the rules, if B = T, 
or T> = T, or the spin of the outgoing line is down, the 
update is rejected. Finally, one has to check whether the 
proposal may lead to the irreducible diagram. For the 
update shown in Fig. 7 the acceptance ratio is given by 
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FIG. 7. Upper panel: An illustration of the Move-P update 
shifting <S along the propagator line. This move changes the 
momentum of the propagator line and flips its spin as well 
as the status of interaction lines between physical and un- 
physical. Lower panel: An illustration of the Move-I update 
shifting S along the interaction line. This move changes the 
auxiliary momentum of the interaction line only. 



Move-I In this self-complementary update one selects 
at random S or T and proposes to shift the selected worm 
along the interaction line to vertex B, see the lower panel 
in Fig. 7. One has to check whether the proposal may 
lead to the irreducible diagram. This update is always 
accepted since it changes only the auxiliary momentum 
variable. 

Commute All possible topologies in a graph of order n 
can be generated by randomly connecting outgoing prop- 
agator lines to incoming ones. Having this observation 
in mind the self-complementary Commute update pro- 
poses to swap destination vertexes for propagator lines 
originating at S or T, see Fig. 8. This proposal is valid 
only if all vertexes have the same space coordinate and 
both propagators have the same spin index. The cru- 
cial advantage of the worm algorithm at this point be- 
comes clear — momentum conservation law is satisfied by 
absorbing the difference k — p into the worm momentum 
Pw ~ * Pw + k — p. This update is accepted with ratio 



ommutc — 



GG 



GG 



(22) 



In the lower panel of Fig. 8 we show a typical diagram 
change produced by Commute. It always changes the 
number of fermionic loops in the graph and, in partic- 
ular, will transform a diagram with a bubble attached 
to the worm vertex into the vertex-correction type dia- 
gram. This explains, to some extent, our design of up- 
dates Create-H and Delete-H. 
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FIG. 8. Upper panel: In the Commute update the diagram 
topology is changed by re-directing propagators originating at 
special vertexes S and T to have their destination vertexes at 
B and A, respectively. Lower panel: typical transformation 
produced by the Commute update. 

Dummy To place the dummy line mark on any of the 

G- or IF-lines we select one of the vertexes at random, 
say vertex A, and make a random decision whether the 
new dummy line should be the interaction line attached 
to A (it has to be of the W type) or the propagator line 
originating from A. The proposal is always accepted if 
the notion of the dummy mark does not change the value 
of the function behind the line (which is assumed to be 
the case here). 

The above set of updates is sufficient for doing the 
simulation. It can always be supplemented by additional 
updates which do not necessarily involve special vertexes 
but reduce the autocorrelation time and lead to more ef- 
ficient sampling of the diagram variables. Moreover, a 
non-trivial check of the detailed balance for debugging 
purposes is only possible if the set of updates is overcom- 
plete. Below we describe several such updates. 

Insert-Remove An idea here is to increase/decrease 
the diagram order by inserting/removing a ladder-type 
structure. More precisely, in Insert we make a random 
choice between S and T to start the construction from 
the special vertex Vi. Next we identify vertex A as the 
destination vertex of the propagator originating from Vi , 
and vertex B as the originating vertex for the propagator 
with the destination vertex V2, which is the other worm 
end, see the upper panel in Fig. 9. H A = Vi the up- 
date is rejected. The proposal is to insert new vertexes 
C (intermediate between V\ and „4) and T> (intermediate 
between B and V2) and to link them with the interaction 
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line of randomly chosen type, either J or W, and momen- 
tum q. The new time variables are generated from the 
probability density t{r). The momenta of new lines and 
the worm are modified as described in Fig. 9 to satisfy 
conservation laws. Finally, if one of the propagators is the 
dummy line, the new dummy line has to have the same 
originating vertex. In Remove we select one of the spe- 
cial vertexes, Vi, at random and verify that the topology 
of lines connecting it to the other special vertex, V2 , as 
well as lines parameters are consistent with Fig. 9 (upper 
panel) . If either C-A or T>- V2 propagator is a dummy line 
the update is rejected. The proposal then is to remove 
vertexes C and T> from the graph and update momenta 
of the lines accordingly. The acceptance ratios are given 
by 



-Rln 



Dy< U W 2 

D v u 9 t(r) 
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(J) 

(W) 



(23) 



vertexes is of a special type) . The proposal is to add new 
vertexes T> (intermediate between B and A) and £ (inter- 
mediate between A and C) and to link them with the W 
line with random momentum q. The new time variables 
are generated from the probability density t(r). The mo- 
menta of new lines are modified as described in the lower 
panel of Fig. 9 to satisfy conservation laws. In Undress 
we select vertex A at random, identify vertexes T>, B, £, 
and C using links along the propagator lines, and verify 
that the topology of lines and their parameters are con- 
sistent with the dressed vertex configuration. If the T>-£ 
line is not of the diagonal W type or one of the propaga- 
tors V-A or £-C is a dummy line, the update is rejected. 
The proposal is to remove vertexes T> and £ from the 
graph. The acceptance ratios are 
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where D„< /D u = \GGGGW/GG\ in Insert and its inverse 
in Remove. 
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FIG. 9. Upper panel: Increasing/decreasing the diagram 
order using a pair of complementary updates Insert and Re- 
move. When the propagators Vi-,4 and B-V2 are linked with 
the new interaction line C-T> carrying momentum q the closed 
loop for momentum conservation goes as Vi-C-U-WVi. The 
same loop is used in the Remove update. Lower panel: Dia- 
gram transformation when vertexes are dressed and undressed 
with interaction lines. 

Dress-Undress One of the easiest updates to increase 
the diagram order within the G 2 M /r -skeleton formulation 
is to dress an existing vertex with interaction line and 
consider the smallest closed loop for transferring momen- 
tum. The Dress update starts from random selection of 
vertex A and identification of vertexes B and C linked to 
it by propagator lines; if B = A the update is rejected 
(we pay no attention in this update whether one of the 
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D„ un (n- 1)*(t)*(t') 



(26) 



where D v > /D u = \GGGGW/GG\ in Dress and its inverse 
in Undress. 

Recolor An easy and efficient way to change the spin 
index of propagator lines is to select a random vertex 
and use it to construct a closed loop by following the 
propagator lines attached to it. If all propagators in the 
loop have the same spin index a it can be changed to 
—a with acceptance ratio unity in the absence of exter- 
nal magnetic field; otherwise, one has to use the ratio 
of products of all propagator lines after and before the 
update. 
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FIG. 10. An illustration of the Move- T update changing the 
imaginary time location of a randomly chosen vertex. One is 
free to change the type of worm from S to T or to place it on 
any of the vertexes in both panels. 

Move- T This self-complementary update is designed 
to sample time variables of the diagram without changing 
its order and topology. The proposal is to select one of 
the vertexes at random, let it be vertex A and update 
its imaginary time variable from r to r' using probability 
density distribution t(r'), see Fig. 10 for two alternatives. 
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The interaction line attached to A cannot be of the J- 
type; whether it is the physical lU-line or unphysical in- 
line does not mater. The acceptance ratio is given by 

_ t{r) j\W(F)/W(F)\\GG/GG\ (generic) 
KMovc ' T ~ t(r')\\F/F\ (bubble)' 

(27) 

D. Diagram sign 

With the dummy line removed, the diagram phase re- 
quired to compute the self-energy and polarization oper- 
ator using Eqs. (12) and (14) is determined by standard 
diagrammatic rules (see also Eq. (10)): 

<p v = arg(/ Une (i/)) + n(n + 1) , (28) 

lines 

where I is the number of fermionic loops (one only needs 
to know whether it is even or odd) . This phase is readily 
recalculated in updates without addressing the whole di- 
agram since I always changes its parity when Create-H, 
Delete-H, and Commute are accepted. 

E. Satisfying the sum rule 

The value of the spin-spin correlation function x( r 
I). ■ ((S* 2 ) 2 ) = 1/4, see Eq. (4), provides an im- 
portant sum rule in the Fourier space 

rV / dq = x(q,n) = V4 , (29) 

which can be used for modifying convergence properties 
of the self-consistent scheme as follows (the integral is 
taken over the Brillouin zone (BZ)). When the maximum 
diagram order is fixed at N the sum rule is violated by 
some amount which vanishes as N — > oo. Since the final 
result is claimed after taking the limit, it is perfectly 
reasonable to impose a condition that the sum rule is 
always satisfied by scaling n by an appropriate factor. 
This is exactly what is done in this article: after solving 
the Dyson Equation we check the value of %(r = 0, r = 
0), adjust the scaling factor for n, and go back to solving 
the Dyson Equation again until the sum rule is satisfied 
with three digit accuracy. 

V. TRIANGULAR LATTICE HEISENBERG 
ANTIFERROMAGNET 

In the diagrammatic formulation there is no concep- 
tual difference in the implementation of the numerical 
scheme for any dimension of space, lattice type, and in- 
teraction range. Thus, sign-problem free systems, e.g. 
the square/cubic lattice Heisenberg antiferromagnet with 



nearest neighbor coupling, can be used for testing pur- 
poses since their properties are known with high degree 
of accuracy (with reliable extrapolation to the thermo- 
dynamic limit) using path-integral and stochastic series 
expansion MC methods. After passing such tests, we 
turn our attention to the triangular-lattice Heisenberg 
antiferromagnet (TLHA) which is a canonical frustrated 
magnetic system with massively degenerate ground state 
in the Ising limit. 

The most important question to answer is whether the 
sign blessing phenomenon indeed takes place, i.e. there 
is a hope for obtaining accurate predictions in the strong 
coupling regime by calculating higher-and-higher order 
diagrams despite factorial growth in the number of con- 
tributing graphs. In Fig. 11 we show comparison between 
the calculated answer for the static uniform susceptibility 

Xu = X (q = 0,m = 0) = f 1 drV X (r,r), (30) 
Jo r 

and the high-temperature expansion results 37,38 at 
T/J = 2. This temperature is low enough to ensure 
that we are in the regime of strong correlations be- 
cause Xu is nearly a factor of two smaller than the 
free spin answer x u = 1/4T. On the other hand, this 
temperature is high enough to be sure that the high- 
temperature series can be described by Pade approxi- 
mants without significant systematic deviations from the 
exact answer 37 ' 38 (at slightly lower temperature the bare 
NLC series start to diverge). We clearly see in Fig. 11 
that the BDMC series converges to the correct result 
with accuracy of about three meaningful digits and there 
is no statistically significant change when more than a 
hundred thousand of 7-th order diagrams are accounted 
for. [We recall that the number of topologically distinct 
diagrams within the G 2 lU-skeleton scheme was calcu- 
lated in Ref. 48; for the eight lowest orders they are 
1, 1, 6, 49, 542, 7278, 113824, 2017881.] The error bar for 
the 7-th order point is significantly increased due to ex- 
ponential growth in computational complexity. The 4-th 
order result can be obtained after several hours of CPU 
time on a single processor. 

Interestingly enough, when temperature is lowered 
down to T/J — 1, which is significantly below the point 
where the bare NLC series start to diverge (see Fig. 13), 
the BDMC series continue to converge (see Fig. 12). 
This underlines the importance of performing simulations 
within the self-consistent skeleton formulation. 

In Fig. 13 we show results of the BDMC simula- 
tion performed at temperatures significantly below the 
mean-field transition temperature. For all points we ob- 
serve extremely good agreement (essentially within our 
error bars) with the Pade approximants used to ex- 
trapolate the high-temperature expansion data to lower 
temperature 38 . Within the current protocol of dealing 
with skeleton diagrams we were not able to go to lower 
temperature due to the development of singularity in the 
response function (and thus effective interaction W at 
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FIG. 11. (Color online) Uniform susceptibility calculated 
within the G 2 W-skeleton expansion as a function of the max- 
imum diagram order retained in the BDMC simulation (black 
dots) for T/J — 2. The result of the high-temperature expan- 
sion (with Pade approximant extrapolation) 38 is shown by red 
square and horizontal line. 
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FIG. 12. (Color online) Uniform susceptibility as a function 
of the maximum diagram order (black dots) for T/J — 1. The 
result of the high-temperature expansion (with Pade approx- 
imant extrapolation) 38 is shown by red square and horizontal 
line. Its error bar is based on the difference between various 
expansion/extrapolation schemes. 



the wave- vector Q = (47r/(3a),0). When the denomina- 
tor 4— J(q)H(q, m) in (7) is close to zero it becomes very 
difficult to control highly non-linear sets of coupled inte- 
gral equations given finite statistical noise on the mea- 
sured quantity Yl(q,m). 

This is clearly seen in Fig. 14 where we show data for 



the staggered susceptibility 0), defined as 

rl/T 



Xs = x(Q,m 



0)=/ 



e^ r x(r,r), (31) 



along with the Curie law and the uniform susceptibility, 
on the double logarithmic scale. 
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FIG. 13. (Color online) Uniform susceptibility as a function 
of temperature (red dots) for the triangular Heisenberg an- 
tiferromagnet calculated within the BDMC approach. NLC 
expansion results 3 ' based on triangles (labeled as 7T and 
8T) and sites (labeled as 12S and 13S) are shown along 
with two different Pade approximant extrapolations of high- 
temperature expansions 38 . 

One of the advantages of our approach is the ability 
to perform calculations of susceptibility at arbitrary mo- 
mentum. In Fig. 15 we show data for x(<l> 0) along the 
r — K — M — r trajectory in the Brillouin zone (BZ). 
Here Y is the center of the BZ, K = Q = (An /(3a), 0), 
and M — (n/a,n/(y/3a)) is the mid-point on the face 
of the hexagonal BZ (see Fig. 15). Results presented 
in Figs. 14 and 15 are new because they are obtained 
for the static (zero Matsubara frequency) susceptibility 
in the thermodynamic limit. It is useful to note that the 
static response is far more difficult to get within the NLC 
method which is suited for calculations of the equal time 
correlation functions, such as, for example, equal time 
spin structure factor. The exception is represented by 
the uniform, or zero momentum, response which is based 
on the total magnetization commuting with the Hamil- 
tonian. In addition, we can afford very high resolution in 
momentum space which is not the case for calculations 
based on clusters of finite size. 

It is clearly seen in Fig. 15 that around T/J = 1 sys- 
tem's response is enhanced along the whole Brillouin zone 
boundary indicative of the frustrated behavior. Only at 
temperatures below T / J = 0.5 it becomes evident that 
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FIG. 14. (Color online) Staggered susceptibility at the wave 
vector Q as a function of temperature (black dots) plotted for 
comparison along with the Curie- Weiss law (blue curve) and 
uniform susceptibility (red dots and line). 
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FIG. 16. (Color online) Modulus of the spin susceptibility 
in real space at T/J = 0.375 on the logarithmic scale. Lat- 
tice points are enumerated according to their distance from 
the origin. Spins on red sites (0, 2, 5, 6, 10) are correlated 
ferromagnetically while spins on blue and black sites are cor- 
related anti-ferromagnetically with the spin at the origin. At 
slightly higher temperature T/J = 0.5 black points (3, 7) 
start to correlate ferromagnetically with the spin So at the 
origin, contrary to the classical ground state pattern depicted 
to the right. 
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FIG. 15. (Color online) Static spin-spin correlation function 
along the characteristic trajectory in the Brillouin zone. 



the system wants to develop correlations commensurate 
with the K point. We confirm previous observation 49 
that even at T/J — 0.375 the spin correlation length, 
which can be estimated from the half-width of the peak 
around K point, is still of the order of lattice constant 
a. This can be checked even more explicitly by looking 
on the static spin correlations in real space. Figure 16 
shows that while the sign structure of short-range spa- 
tial correlations is consistent with the three-sublattice 
120° state, the magnitude of the correlations becomes 



exponentially small on the scale of a few lattice periods. 
Closer look also reveals dramatic suppression of correla- 
tions between site (say, sublattice A) and sites 3 and 7 
both of which would belong to sublattice C in the per- 
fectly ordered 120° classical state. Moreover, at slightly 
higher temperature T/J = 0.5 the sign of correlations 
on sites 3 and 7 changes sign and turns ferromagnetic, 
similar to A-sublattice spins, though with much smaller 
amplitude. This temperature induced reversal of cor- 
relations is a remarkable effect specific for a frustrated 
system. 

It has been noted some time ago that short wavelength 
spin excitations contribute significantly to the finite tem- 
perature properties of triangular lattice antiferromagnet 
at not too low temperature 50 . This has to do with sub- 
stantial phase space volume these excitations occupy as 
well as with their relatively weak dispersion 51 . It is con- 
ceivable that particularly weak correlations between sub- 
lattices A and C noted above have to do with these ex- 
citations as well. All of these features can be extracted 
from the retarded spin susceptibility xi.1i - 1 ) calculation 
of which requires analytic continuation of our Matsubara- 
frequency susceptibility xili^m) to the real frequency. 
We plan to address this important issue in the near fu- 
ture. 



VI. CONCLUSIONS 

This paper describes novel approach to frustrated spin 
systems. Obtained numerical results for the spin- 1/2 tri- 
angular lattice Heisenberg model, Section V, show the 
power and competitiveness of our approach in compari- 
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son with other well established numeric techniques. 

Future work has to address the issue of performing 
simulations at lower temperature in the regime charac- 
terized by the large correlation length. Technically, this 
translates into being close to zero in the denominator 
of Eq. (7) for some momentum values. Progress in this 
direction should allow us to better describe the cross- 
over/transition from the cooperative paramagnet to the 
long-range ordered (albeit frustrated) state. 

Yet perhaps the most promising line of attack has to 
do with applying our technique to the geometrically frus- 
trated models that do not support magnetically ordered 
state at all. In two dimensions this singles out quantum 
kagome lattice antiferromagnet which have recently be- 
ing shown to realize a long-sought Z2 spin liquid state 27 . 
This task will require extension of our approach to sys- 
tems with several (three in this case) spins in a unit cell. 
The added matrix complexity does not represent any fun- 
damental difficulty. 

Moving one dimension higher brings one to the most 
frustrated antiferromagnet in the world - spin-1/2 py- 
rochlore antiferromagnet 52 . Pyrochlore Ising-like 
model realizes beautiful quantum spin-ice physics 55 while 
the fate of spin-1/2 Heisenberg model is an open ques- 



tion. It is widely believed that 'cooperative paramagnet' 
region is most extended in this three-dimensional system. 
Unlike many lower-dimensional frustrated system, spin- 
1/2 pyrochlore is essentially not accessible by quantum 
Monte-Carlo technique due to large unit cell (4 spins). 
We believe that our Diagrammatic MC approach is there- 
fore uniquely suited for studying the finite-temperature 
dynamics of this outstanding frustrated magnet. 

Recently we have generalized the Popov-Fedotov trick 
to a universal technique of fermionization which leads to 
a well-defined standard diagrammatic technique for arbi- 
trary lattice spin, boson, and fermion system with con- 
straints on the on-site Fock states 39 . This development 
creates a broader context for the present work: success- 
ful implementation of the BDMC method for models of 
quantum magnetism may lead to the universal numer- 
ical tool for arbitrary strongly correlated lattice models 
within the fermionization framework when diagrammatic 
expansion do not involve large parameters. 

We thank M. Rigol for communicating us data ob- 
tained within the NLC method. This work was sup- 
ported by the National Science Foundation under grants 
PHY-1005543 (S.K., N.P., B.S., and C.N.V.) and DMR- 
1206774 (O.A.S.), and by a grant from the Army Re- 
search Office with funding from the DARPA. 
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